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SUMMARY 


The theory of generalized matrix inversion is utilized in formulating a 
recursive algorithm for least squares parameter estimation. This algorithm 
allows the parameter estimation to Begin after the first "observation" has Been 
made and affords a means of computing the nth parameter state from the (n-l)st 
"parameter state" and the nth "observation." The problems associated with sin- 
gular matrices encountered in iterative least squares procedures do not affect 
the algorithm. 


INTRODUCTION 


In the theory of linear least squares parameter estimation, the matrix 
equation Ax = b is encountered, where A is an n X k matrix, x is a 
k X 1 "parameter" state vector, and b is a n X 1 "observation" vector. 
There is usually no vector x that will satisfy this matrix equation so that 
in some sense a "best solution" is to be found. In the least squares theory 
the "best solution" is defined to be that k x 1 vector x such that 


f(x) = (Ax - b) T (Ax - b) is minimum. It is well known that if x 
f(x), then x must satisfy the normal equations A Ax = A b. If 


minimizes 

T 

A A is 


nonsingular and f(x) has a minimum then, indeed, that minimum is attained at 

^ / rp \-l ip 

x = \A Ay A b. It will later be shown that, in any case, the matrix equation 
T T 

A Ax = A b always has some solution and, in fact, may have infinitely many 

T 

solutions. In case A A is singular, a particular solution that has physical 
meaning can be chosen. 


Nonlinear parameter estimation problems are usually handled by linear ap- 
proximations of the actual parameter state in a neighborhood of a nominal pa- 
rameter state. The resulting equations are of the same general form Ax = b; 
however, in this case x denotes the deviation from the nominal state, and b 
denotes the deviation in "observed" and "computed" values. The solution to this 
equation requires iterative procedures and involves problems of singular matri- 
ces. The recursive algorithm to be developed will not be affected by these sin- 
gularity problems. 
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The recursive algorithm differs from a similar algorithm developed by 
Gainer (ref. l) in that no matrix inversion is required as a result of waiting 
for a sufficient number of observations to accumulate. The estimation proced- 
ure may begin after the first observation is made. The parameter state can be 
estimated at a fixed epoch in time (in dynamic systems) and, hence, is useful 
in trajectory calculations as well as in guidance procedures. There is no need 
for matrix inversion (and associated storage requirements ) in onboard calcula- 
tions . 


SYMBOLS 


Capital letters 
Lower case letters 
Greek letters 



diag a , 
Z 
0 



a 


n 


matrices 

vectors (unless otherwise stated) 
scalars 

transpose of A 
inverse of A 
generalized inverse of A 
diagonal matrix 
zero matrix 
zero vector 
euclidean norm 


THE GENERALI ZED INVERSE 


A. Bjerhammar (ref. 2 ), E. H. Moore (ref. 3)> and R. Penrose (ref. 4) inde' 
pendently generalized the concept of matrix inversion. The generalized inverse 
of a singular and nonsquare matrix possesses properties which make it a central 
concept in matrix theory. Only real matrices are considered in the definitions 
and theorems to follow. However, for complex matrices, the definitions and 
theorems are identical if the word "transpose" is replaced with "conjugate 
transpose." The following fundamental theorem due to R. Penrose (ref. 4) is 
stated without proof. 
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THEOREM I. The four equations 


AXA 


= A 


XAK = X 
(XA.) T = XA 
(AX) T = AX 

have a unique solution X for each real matrix A. 


(1) 

( 2 ) 
(3) 

oo 


The solution X in THEOREM I is denoted as X = A and called the gener- 
alized inverse of A. It is easy to see that the defining equations for A + 
imply that AA + and A + A are, respectively, orthogonal projection operators 

on the range spaces of A and A + . For the sake of completeness the next the- 
orems vill give some veil knovn properties of the generalized inverse. 

THEOREM II. Let A be an arbitrary real matrix. Then, for scaler \ ± 0 
and unitary U and V 


a + (a + ) t a t -a + -a t (a + )V 


a + aa t = a t = a t aa + 


(a + ) + > A 

(a T ) + - (a + ) t 


A + = A" 1 


for nons ingular A 


( 5 ) 

(6) 

(7) 

(8) 

(9) 


(AA) + = ^A + 

(a t a) + = a + (a + ) t 

(UAV) + = V^aV -1 

A = > A. and aTa . = Z 

L i 1 3 

A^A. = Z for i / j 

3 i 

imply 



( 10 ) 

( 11 ) 

( 12 ) 

(13) 
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If A is normal (i.e. A^A = AA^) 
then, A + A = AA + and (A n ) + = (a + ) n 


T 4" 

A, A A, A and A A all have rank equal 
to trace 


A + A 


- (a t a) 


+.T 


(14) 

(15) 

(16) 


Note that equation (l6) reduces the problem of computing A to that of 

T 

computing the generalized inverse of a symmetric matrix A A. Moreover, such 
a matrix can always he diagonalized by a unitary transformation, that is, 


D = u(a T a)v = diag . 
Now equation (12) implies that 


* n ) 


(a t a) + - 


VD U = V diag - , 


-i}> 


diag - , 


It is tacitly assumed that if a. = 0, the corresponding term in 
h ' i \ 1 

. , — is zero. It is not usually an easy task to determine the 
v 1 &n / 

unitary transformations U and V. Methods for computing the generalized in- 
verse have been given by various authors (refs. 2, 5, 6, 7 > and 8). 

The following is a theorem of major importance, characterizing all solu- 
tions of the matrix equations AXB = C which have some solution X. 

THEOREM III. For the matrix equation AXB = C to have a solution, a necessary 
and sufficient condition is 

AA + CB + B = C 

in which case, the general solution is 


X = A + CB + + Y - A + AYBB + 


where Y is arbitrary to within the limits of being consistent with dimension 
in the indicated multiplications (ref. 4). 

Proof: Suppose X satisfies AXB = C. Then, 

C = AXB = AA + AXBB + B = AA + CB + B 
+ + + + 

Conversely, if C = AA CB B, then A CB is a particular solution. Clearly, 
for the general solution AXB = 0 must be solved. Any expression of the form 
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X = Y - A + AYEB + 

is such a solution. Moreover, if AXB = 0, then, 

X = X - A + AXBB + 

Hh "4* 4* 

The only properties required of A and B in the theorem are AA A = A and 
BB + B = B. 

Corollary A. The general solution to the vector equation 

Px = c 


is 

X = P + c + (i - P + p)y 

where y is arbitrary, provided a solution exists. 

Corollary B. A necessary and sufficient condition for the equations 

AX = C 


to have a common solution is that each have a solution and AD = CB (ref. 4). 

Proof: If AX = C and XB = D have a common solution, then clearly each has 

a solution, and 

AXB = CB 
AXB = AD 


so that 


CB = AD 


In order to obtain the sufficiency, it is assumed that 

X = A + C + DB + - A + ADB + 

which is a solution if AD = CB, AA + C = C, and DB + B = D. 
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THEOREM TV. The terms 


and 


A + A, AA + , I-A + A, and I - AA + 


are symmetric idempotents 


( 17 ) 


H is a symmetric idempotent which implies that 


H = H 


(18) 


THE RECURSIVE ALGORITHM 


An easy consequence of equation (6) is that the matrix equation 

T T + 

A Ax = A b always has a solution. Indeed, x = A h is a solution. In fact, 

T 

if A A is nonsingular, this solution is unique and is the usual least squares 

solution. Moreover, according to corollary A, x = A + b + (i - A + A)y, for 
arbitrary y, gives all solutions. The triangle inequality also implies that 

x = A + b is that solution with minimum norm. When x is a parameter state 

deviation vector (i.e. the nonlinear estimation case), x = A + b yields the 
estimate of the parameter state deviation that has minimum norm. Physically 
speaking, this means that among all possible estimates of the actual parameter 
state, the estimate of the actual parameter state obtained by correcting the 
nominal parameter state by x is as close as possible (in norm) to nominal 
parameter state. 

T T 

The possibility of A Ax = A b having infinitely many solutions is 

T 

clearly equivalent to the singularity of the matrix A A. In this case an 
iterative procedure would halt. However, the adroit choice of the physically 
meaningful minimal norm solution would allow the iteration to continue. 


From a computational viewpoint, it is desirable to have a recursive 
algorithm for the state estimation. The following theorem will give rise to 
the recursive computation of the parameter state deviation. 


THEOREM V. 


Let 



Hh 

a 


and 



( 19 ) 


( 20 ) 
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I 


f 

l 


4- XT 

Proof: It will be shewn that a - a satisfies equations (l) to ( 4 ). 

a 



4- T 

Greville (ref* 9 ) noted that a is a constant multiple of a but does not 
give this explicit form. 


As a simple example, consider the vector equation 

(i,D (*) - (2) 

This equation has a solution a = 1, b = 1. In fact, the equation has infinitely 
many solutions so that according to Corollary A the general solution is given by 



- ( 1 , 1 ) + ( 2 ) + 



- ( 1 , 1 ) + ( 1 , 1 ) 


y 


vhere y is an arbitrary 2X1 
components of y and note that 


vector. Let c 



and d denote the arbitrary 
to obtain 



Since c and d are arbitrary, c-d = f is arbitrary, and the following 
equation results: 
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These coordinates (see fig. l) describe a straight line whose equation is 
h = 2 - a. Hence, any vector whose coordinates satisfy b = 2 - a (e.g., R, 
whose end point lies on the line) is a solution of the given vector equation. 

Note that the solution with minimum norm is R Q = (l,l) . In a practical sense, 

the origin might be considered to be a nominal set of initial conditions, a 
given nominal parameter state, and so forth. Figure 1 gives the geometrical 
significance of the minimum norm solution. 
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THEOREM VI. The following theorem is 
due to G-reville (ref. 9 )- Let F he 
any m X p matrix, f denote the pth 


column of F, and F ^ represent the 

suhmatrix consisting of the first p - 1 
columns of F, that is. 


F = 




then, letting 


I 




f 

P 


s 

P 


the following equation results: 



where 


c 

P 



f 

P 


d 

P 


(V 



T \ -1 T + 
c c c F 
PP P P-1 



* e 



= e 


THEOREM VII. Let A be any n X k matrix, a denote the nth row of A, 
and A^ ^ represent the submatrix consisting of the first n-1 rows of A, 
that is, 
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A = 



Then, letting 


a (i - A + .A ^ = h 
n \ n-1 n-1/ n 

the follcrwing equation results 

A+ ‘( A n-1 - *zfln : *n) 

■where 


q = a A n 
ha n n-1 


n 


p n = ( 


.T\" 1 .+ T 


1 + VJ VA 


Proof: If 


'A , 

a =i .?:! 

Et 

n 


then according to equation (8) of THEOREM II 


A + = 


.T T 

A _ : a 
n-1 • n 



Clearly 



n-1 



f 

n 



has the form described in THEOREM VI vith n = p, m = k, 

and s = (i - A T A 14- A a T = (i - A + A la T . 
n y n-1 n-1/ n l n-1 n-1/ n 


Hence , 
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where 


c 


n 



d 

n 



T 

c c 
n 


O' 


T.T+ 
c A . 
n n-1 





e ) 

6 ) 


Noting that 
equations result: 
T 


A + = ( A + , - a T c T : a T ) 

\ n-1 n n • n/ 

«r - 


= h and s 
n n 


= 0 if and only if h = 0, the following 


n 


c = a A _ 
n n n-1 


d T 

n 


, T V 1 
+ c c , 
n n/ 


.+ 

A c 
n-1 n 



e ) 



9 ) 


so that defining 



T 

and q^ = c^ completes the proof of the theorem. 


Using the results and notation of THEOREM VII, 
of Ax = h can he realized as a recursive process, 
noted that the least squares solution after all n 
given hy 


the least squares solution 
To this end it need only he 
observations are made Is 


A 



= A 


b 


+ 

so that writing A h 


in partitioned form the following equation is obtained: 


x 

n 




where 

1X1 


h _ is a column vector of the 
n-1 

vector corresponding to the nth 


first n-1 observations and 
observation. Multiplication 


b is a 
n 

yields 


\ 

\ 
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S n ■ VlVl - V^n-l + p n b n 

- S n-1 - VtAi-I + p n b n 


Nov since 


<n 


a A , the last expression may he written as 
n n-1* 


x = 


= x - p a x _+ph 
n-1 n n n-1 n n 

= (i - p a ^x n +ph 
\ n/ n-1 ^n n 

Note that this recursive algorithm- involves at most the generalized inver- 
sion of the row vector a ( I - A + n A . ) . All other quantities are either 

n\ n-1 n-1/ ^ 

known from the (n-l)st state, or they are simple functions of the nth obser- 
vation. The generalized inversion of the row vector - 'kn-l^'n-l) 

trivial in light of THEOREM V. 

THE COVARIANCE OF THE ESTIMATE 

Consider the vector equation Ax - b = e where e is an n X 1 error 
vector. In the absence of weighting it is usually assumed that E(e) = 9 and 

E(ee^) = I, where E denotes the expected value operator. The covariance 
matrix, C^x n ,x^, of the nth estimate is given by 


(tA) * E (% * x ) ( - x ) 


T 


The equation Ax - b = e may have infinitely many solutions; however, the 
minimal norm solution x = A + (h + e) is chosen. With this solution and the 
least squares estimate x^ = A^h, the following equation results: 


C (V*n) = E R A+t " x )( A+b - x ) 

= E [(A + e)(A + e) T ] 

= E (A + ee T A +T ) 

= A + E (ee T ) A +T 


T 


= A + IA +T = A + A +T 
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The recursive computation of cfx^x^J from C^x^ ^ and the nth 

observation is achieved in the following way: 

°(V S n) ■ A+A+T ■ ( A n-1 - Vn I p n) ( A n-1 ' Vn : p n) 

* ( A n-1 - 5 A)(Vl ' p n%) * p n p n T 

,+ +T .+ T T ,+T . T T 

■ A „-lVl - Vl Vn - p nVn-l + p nWn 


Since 


\ =># 


+ 

' -Li- 

n n-1 


then. 


C(x ,x ) = C(x , ,x \ - C(x n ,x -^a^p 
\ n’ rij V n-1* n-ly \ n-1* n-V n^n 

/■n /^* /\ \ „ /a /x \ T T 

-pa Cfx _,x _ ) + p a Cfx n ,x _ )a p 

^n n \ n-1' n-1/ ^n n V n-1* n-1/ n^n 


WEIGHTED OBSERVATIONS 


T — 1 

To minimize (Ax - b) VT (Ax - b), that is, 'weighted least squares, note that 
W is usually a positive definite symmetric covariance matrix and hence there 

T -1 ~ ~ 

exists a matrix Q such that Q Q = W . For A = QA and b = Qb the theo- 
retical results are the same. 


CONCLUDING REMARKS 

The recursive least squares equation, "by its very nature, does not Involve 
singularity problems. In nonlinear parameter estimation procedures the minimal 

norm solution of the parameter state deviation is always given by A + b. The 

T 

minimal norm solution Is the only solution when A A is nonsingular. However, 

T 

in iterative procedures the singularity of A A would not bring the iteration 
to a halt. In any case, the state estimation may begin after the first obser- 
vation is made. An additional computational feature of the algorithm is a small 
machine storage requirement. Note that the computation of x^ only requires 

storage ffrom computation of x^ ) of the k X k matrices A* ^A^ and 
+ +T zs, ^ n n 

A _ A n and the estimate x _ . Such a feature would be important, for in- 
n-1 n-1 n-1 ^ 9 

stance, in onboard spacecraft computations. 
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